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"^o" ■ Abstract 
O 

In this paper we propose a microscopic model to study the polymerization of microtubules (MTs). Starting 

(N ; 

, from fundamental reactions during MT's assembly and disassembly processes, we systematically derive a non- 

O . 

' linear system of equations that determines the dynamics of microtubules in 3D. We found that the dynamics 
• of a MT is mathematically expressed via a cubic-quintic nonlinear Schrodinger (NLS) equation. Interestingly 
the generic 3D solution of the NLS equation exhibits linear growing and shortening in time as well as temporal 
fluctuations about a mean value which are qualitatively similar to the dynamic instability of MTs observed 
experimentally. By solving equations numerically, we have found spatio-temporal patterns consistent with 
experimental observations. 
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I. INTRODUCTION 



More than 25 years ago, Del Giudice et al. [1] argued that a quantum field theory approach to 
the collective behaviors of biological systems is not only applicable but the most adequate as it leads 
naturally to nonlinear, emergent behavior which is characteristic of biological organization. They 
followed a line of reasoning championed by Davydov [2, 3] and Frohlich [4] who emphasized integration 
of both conservative and dissipative mechanisms in biological matter leading to the emergence of spatio- 
temporal coherence with various specific manifestations such as almost lossless energy transport and 
long-range coordination. 

Microtubules (MTs) are long protein polymers present in almost all living cells. They participate 
in a host of cellular functions with great specificity and spatio-temporal organization. Microtubules 
are assembled by tubulin polymerization into a helical lattice forming a cylinder which is rigid and 
straight by biological standards. These protein polymers, typically several microns long, participate 
in fundamental cellular processes such as locomotion, morphogenesis, and reproduction [5]. It is also 
suggested that MTs are responsible for transferring mechanical energy across the cell, with little or no 
dissipation. 

Both in vivo and in vitro observations confirmed that an individual microtubule switches stochas- 
tically between assembly and disassembly states making MTs highly dynamic structures [6]. This 
property of MTs is referred to as dynamic instability and it is a nonequilibrium process. It is gener- 
ally believed that the instability starts from the hydrolysis of guanosine triphosphate (GTP) bound to 
tubulin converting it to guanosine diphosphate (GDP). This reaction is exothermic and releases ~ 8kT 
energy per reaction [7], i.e. approximately 0.22 eV per molecule [8] at room temperature. Here k is the 
Boltzmann constant and T is the absolute temperature. Since GDP-bound tubulin favors dissociation, 
a MT enters the depolymerization phase as the advancing hydrolysis reaches the growing end of a MT. 
This dynamic phase transition is called a catastrophe. As a result, MTs rapidly disassemble releasing 
GDP-tubulin in the solution where reverse hydrolysis takes place followed by a re-polymerization phase 
of MTs called a rescue. Therefore, MTs constantly fluctuate between growth and shrinkage phases. 

Several theoretical models have been proposed for a macroscopic description of these processes 
using nonlinear classical equations [9-20]. These phenomenological models provide good agreement 
with experiment but generally consider a MT as a one-dimensional mathematical object that switches 
stochastically between growth and shrinkage states. Even two-dimensional [9, 10, 19] and three- 
dimensional [20] considerations for MT polymerization have started from a given 2D (planar) or 3D 
(cylindrical) mathematical configuration that can grow or shrink randomly. 

In this paper, however, we based our model on fundamental biochemical reactions that are occurring 
during microtubule's assembly and disassembly processes. This enables us to derive a nonlinear system 



particular, our model explains the continuum symmetry breaking of an isotropic pool of tubulin dimers 
that leads to formation of experimentally observed 3D structures such as ring-shaped or filaments [21- 
23]. We believe that this treatment is necessary to address the fundamental issues about the observed 
dynamical behavior of MTs. As stated by Del Giudice et al. [1] : "Systems with collective modes 
are naturally described by field theories. Furthermore, quantum theory has proven to be the only 
successful tool for describing atoms, molecules and their interactions." 

II. THE MODEL 

Consider an individual MT in a free tubulin solution containing a large number of GTP-tubulin, 
GDP-tubulin and a pool of free GTP molecules. In this solution several processes take place (as well 
as their reverse reactions): 

(i) GTP hydrolysis: 

GTP — ► GDP + Ai, 

(ii) conversion of tubulin GDP from tubulin GTP: 

Tqtp — ► T G dp + A 2 , 

(iii) growth of a MT: 

A3 + MTjv_i + Tgtp — > MT N , 

(iv) shrinkage of a MT: 

Tgdp + A 4 . 

Experimental studies determined the free energy values for these reactions as: Ai ~ 220 meV, A 2 — 160 
meV and A 3 ~ A 4 ~ 13 x 40 = 420 meV, respectively [24]. They are clearly above the thermal energy 
at room temperature (kT ~ 26 meV) and within a quantum mechanical energy range that corresponds 
to the creation of one or a few chemical bonds. Hence we may consider each chemical reaction as a 
quantum mechanical process [25]. As a result, an individual microtubule with length L can be viewed 
as consisting of N tubulin layers defining its quantum state \N). A tubulin layer consists of at least 
one tubulin dimer and at most 13 tubulin dimers as observed in the MT's structure. The state can be 
raised/lowered by a creation/annihilation operator (i.e. polymerization/depolymerization process) to 
the I N+ 1)/| N— 1) state. The corresponding MT is then longer/shorter by one tubulin layer compared 
to the original one. 

Further simplification will be achieved by combining the above processes into two fundamental 
reactions: 

(I) growth of a MT by one dimer length by adding of one tubulin layer in an endothermic process: 



(II) shrinkage of a MT by one dimer length due to the removal of one layer of T G dp dimers in an 
exothermic process: 

MT N — , MTjv-i + Tqdp + A, (2) 

where A is the energy of the reaction. In order to derive a quantum mechanical description of mecha- 
nisms (1) and (2), we introduce quantum states of a MT, tubulin and the heat bath, respectively: 

• | TV) is the state of a microtubule with N dimers (containing both GTP and GDP tubulins). 

• \N T ) is the state of a tubulin dimer, T GTP or T GDP . 

• | TV) is the GTP hydrolysis energy state. 

Then, the relevant second quantization operators would be: 

<f = \N + 1)(N\, a = \N-l)(N\, 

&t = \N T + 1)(N T \, b = \N T -l)(N T \, (3) 
= \N + 1)(N\, d=\N-l)(N\. 

Here b/tf and d/(fi are annihilation/creation operators of tubulin and energy quanta, respectively. 
The operators a/a^ are lowering/raising the number of tubulin layers in a MT. The creation and 
annihilation operators obey the Bose-Einstein commutation relations 

[Qk, qL\ = <*km, and [q[, q£j = = [q k , q m ], (4) 

where [A, B] = AB — BA is the Dirac commutator or q = a, 6, and d . Following [26], one can express 
the above processes using creation and annihilation operators (3): 

a)b d : A + MT A r_ 1 + T GTP — ► MT^ (5) 
d ] b ] a : MT N — ► MT W _! + T GDP + A (6) 

Operators (5) and (6) describe a MT's growth and shrinkage by one layer, respectively. The poly- 
merization or depolymerization process may happen repeatedly before reversing the process which can 
be captured by constructing product operators, i.e. (a^b d) m and {$ a) n , where m and n are the 
number of growing or shrinking events in a sequence, respectively. Based on the mechanisms in (5) 
and (6), the Hamiltonian for interacting microtubules with T G tp/T G dp tubulins can be written as 

H = ^ ^kflkflk + ^2 ^m&L&m + ^2 ^ d l dl 
k m 1 

oo 



where u, w, a and Y are constants that can be related to the energy of the fundamental processes 
[27]. A growing/shrinking MT may change its state quickly or after several steps to a depolymeriz- 
ing/polymerizing state and then may change back to a polymerizing/depolymerizing state. Experi- 
mentally at a mesoscopic level the transition from the growing to the shrinking phase is quantified 
by the catastrophe rate / cat and the transition from the shrinking to the growing phase is expressed 
by the rescue rate / rcs in which / res < / ca t- In the Hamiltonian (7), these transitions correspond to a 
combination of creation and annihilation operators as the n th power of the reaction in (5) and (6): 

c k„m„ = (ati 6 ^ d h )(a{ 2 b m2 d h ) . . . (a{b mn d ln ). (8) 

Here k n = {k 1; k 2 , . . . , k n } is a collection of indices and J2^ n = X^k 2 ■ ■ ■ Sk„- We no ^ e ^ na ^ the 
momentum conservation for the last two terms in the Hamiltonian (7) requires that 1„ = J2i=i k« _ 
Y^i=i m i ~ Y17=i lj- Therefore, the first n — 1 of 1 will be free and summed in the Hamiltonian (7). 

III. THE DYNAMICAL EQUATIONS 

The Heisenberg equation of motion for a space- and time-dependent operator q(r,t) reads: 
ihd t q(r,t) = —[H,q(r,t)], where H is the Hamiltonian. A system of coupled equations that describes 
the quantum dynamics of a MT can be derived from the Heisenberg equation. However, since MTs are 
overall classical objects (although some of their degrees of freedom may behave as quantum observ- 
ables), we need to ensemble average over all possible states to obtain effective dynamical equations. 
Fourier transforming a v , b v and d v operators over all states, we find 

^(r, t) = IT 1/2 exp(-i»7 • r)a„(t), (9) 
v 

x (r, t) = Q- 1 / 2 exp(-i»7 • r)6„(f), (10) 
v 

0(r,O = n- 1 ' 2 J2 e M-iV ■ r)d v (t), (11) 
v 

where Q is the volume over which the members of the plane wave basis are normalized [28, 29]. Here 
ip(r,t), x( r ,^), and (f>(r,t) are the corresponding field operators for the quantum operators a v , b v 
and d v , respectively. The derivation of the equation of motion for the field operators is lengthy but 
straightforward and given in detail in [30]. The final form of the equations of motion is found to be 

idtip + iv • = -6V 2 ^ + Vif), (12) 
id tX = -eV 2 X + U X , (13) 
\/(|^|,|x|) = a + c|x| 4 |^| 2 -rf| X | 6 |^| 4 , 
Um) = f-h\*P\ 2 , 



TABLE I: The parameters available in the literature. 



Parameter 



Simulation Coeff. 



Exp. Value 



Reference 



MT growth rate 
MT shortening rate 
MT catastrophe frequency 
MT diffusion constant 
Tubulin diffusion constant 



Real(c) 
Real(d) 



0.50 - 19.7 (/xm/min) [31] 

4.1 - 34.9 (/xm/min) [31] 

0.12 - 3.636 (/min) [31] 

2.6 - 30.3 (/xm 2 /min) [32] 

300 - 480 (/an 2 /min) [33] 



For simplicity we assumed that the energy in the system is distributed uniformly during the course 
of experiment. As a result, = and V0 = 0. Here parameters a, b, e and / are real but c, d and 
h are complex. Table 1 lists experimental values for some of these parameters. Eq. (12) represents 
the nonlinear cubic-quintic Schrodinger (NLS) equation with a complex potential that has been ex- 
tensively studied in connection with topics such as pattern formation, nonlinear optics, Bose-Einstein 
condensation, superfluidity and superconductivity, etc. [34]. A general solution of the NLS equation 
can be cast in the form of 

1>(T,t)=R(T,t)exp[iS(T,t)], (14) 

which involves topological defects (point in 2D and line in 3D). In 3D space these defects represent ID 
strings or vortex filaments [34] . Furthermore, it is shown that the symmetry group of the cubic-quintic 
NLS equation is an extended Galilei group that includes translational and rotational symmetries as 
well as proper Galilei boosts and total mass conservation [35] . Adding any constrains such as boundary 
conditions to the equation, however, will cause a symmetry breaking in the system. As an example, in 
a cylindrical coordinate system, due to the rotational symmetry breaking, the general solution reduces 
to a stationary solution that represents a straight vortex filament with a twist: 

■0(r, 9, z, t) = R(r) exp[i(ut + n6 + w{r) + k z z)], (15) 

where uj is the spiral frequency, R(r) is the amplitude, w(r) is the spiral phase function and integer 
n is the winding number of the vortex [35, 36]. The axial wave number k z characterizes the vortex's 
twist. In the case of the NLS equation, a family of vortices that move with a constant velocity is also 
a solution [36]. 

Further symmetry breaking would lead to different 3D structures such as double-wall , ring-shaped, 
sheet-like, C-shaped and S-shaped ribbons, and hoop structures as seen during tubulin polymerization 
experiments [21-23]. 



IV. NUMERICAL RESULTS 



Equations (12) and (13) are solved numerically with a no-flux boundary condition. As an initial 
condition we chose a straight vortex filament perturbed by small noise (eg. thermal or environmental 
noise). In Fig. 1 we compare the observed data on the MT length as a function of time with our 
simulation results. The length of a vortex is defined as [36] : 



where 6(x) is the step function and ip is a constant. In Figs. 1 and 2 we compare the observed data 
on the MT length as a function of time with our simulation results. Experimental panels in Figs. 1 
and 2 represent the experimental data published by Rezania et al. [37]. Simulation panels show the 
numerical results of the normalized vortex length as a function of time for the given set of parameters. 

To provide a simple yet accurate and powerful comparison between experimental and simulation 
results, we graph recursive maps for the data points. The advantage of the recursive maps is the 
introduction of regularity into the data sets that allows for a better choice of adjustable parameters 
due to noise reduction inherent in the separation of data into subsets corresponding to independent 
processes. In spite of being very simple, recursive maps of assembly and disassembly processes of 
individual MTs can successfully reproduce many of the key characteristic features. Consider first the 
following stochastic map as the simplest case that illustrates the approach taken: 



where £(t n ) is the length of a microtubule after n time steps, t n . The parameter r is chosen to be a 
random number with the following two possibilities: 



In terms of the MT polymerization process, p is the probability that a given event will result in 
assembly while 1 — p is the probability of a complete catastrophe of the MT structure. The above 
simplified model, therefore, is governed by only two adjustable parameters: (a) the probability of 
complete catastrophe 1 — p which is constant and independent of the length or time elapsed and 
(b) the rate of polymerization which is proportional to the length increment 5 over the unit of time 
chosen in the simulation. Thus, the coefficient 5 divided by the time step At (= t n+ i — t n ) gives the 
average growth velocity of an individual MT. Such information can be used to fine tune the simulation 
parameters. We note that the slope of the line in the simulation panels in Fig. 1 can be adjusted by 
varying the real parts of parameters c and d. The frequency of catastrophe events can also be changed 




(16) 



£(t n + l)=r[£(t n ) + 5], 



(17) 
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FIG. 1: Length of a distinct microtubule as a function of time. The top-left panel represent experimental 
data published in [37]. provided by O. Azarenko and M.A Jordan from the University of California, Santa 
Barbara. The top-right panel is the simulation result with the set of parameters (a = 1, b = 10, c = 10 + i, d = 
20 + i,e = 300, / = 1 and h = + The bottom-left panel shows recursive maps for both experimental and 
simulation results. The bottom-right panel represents the power spectrum of the experimental (solid curve) 
and simulated (dashed curve) data, respectively. The curves are plotted at different offsets for clarity. No 
particular frequency of oscillation can be seen from the power spectrums. However, both power spectrums 
show a very similar broad distribution that more or less decays with frequency as an inverse power-law with 
slope ~ 1.0. The best fit inverse power-law is shown by a dotted line. 

map for both the experimental data and the simulation results. Based on the recursive maps, the key 
characteristics of the experimental and simulated results that were obtained independently are quite 
similar. This represents Eqs. (12) and truly describes the dynamics of MTs' polymerization. 

To provide a more solid comparison, a spectral analysis is also carried on both experimental and 
simulated data. As discussed by Odde et al. [38], the power spectrum analysis is a more general way 
to characterize the microtubule assembly/disassembly dynamics without assuming any model a priori. 
The power spectrum panels in Figs. 1 and 2 represent the spectral power of the experimental and 
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FIG. 2: Same as Fig 1. but with the set of parameters (a = 1, b = 30, c = 10 + 10*, d = 20 + 10*, e = 300, / = 1 
and h = — .1 + i). In the power spectrum panel, the inverse power-law has a slope of ~ 1.2. 

curve) and simulated (dashed curve) spectrums. We note that the curves are plotted at different offsets 
for visual clarity only. As expected, no particular frequency of oscillation can be found from the power 
spectrums. However, both power spectrums demonstrate a very similar broad distribution that more 
or less decays with frequency as an inverse power-law with slope ~ 1.0 in Fig. 1 and ~ 1.2 in Fig. 2, 
respectively. The best fit inverse power-law is shown by a dotted line in both panels. 

Although this is not unique to the model presented here, our results show no attenuation states 
during MT polymerization (Figs. 1 and 2). The MT length undergoes small fluctuations all the 
time. This can be understood by noting that our model is based on the cyclic polymerization and 
depolymerization of tubulin dimers. Behavior consistent with this result has recently been observed 
by Schek et al. [39] who studied the microtubule assembly dynamics at higher spatial (~ 1-5 nm) and 
temporal (~ 5 kHz) resolutions. They found that even in the growth phase, a MT undergoes shortening 



V. DISCUSSION 



The basic structural unit of a MT is the tubulin dimer. Each dimer exists in a quantum mechanical 
state characterized by several variables, especially GTP/GDP. Each microstate of a tubulin dimer 
is sensitive to the states of its neighbors. Tubulin dimers have both discrete degrees of freedom 
(distribution of charge) and continuous degrees of freedom (orientation). A model that focuses on 
the discrete variables will be an array of coupled binary switches [40], while a model that focuses 
on the continuous ones could be an array of coupled oscillators [41, 42]. Here we have focused on 
tubulin binding and GTP hydrolysis as the key biochemical processes determining the states of MTs 
which are most easily accessible to experimental determination. We have shown how a quantum 
mechanical description of the energy binding reactions taking place during MT polymerization can 
lead to nonlinear field dynamics with very rich behavior that includes both localized energy transfer 
and oscillatory solutions. 

Based on the chemical binding reactions occurring during MT polymerization, a quantum mechan- 
ical Hamiltonian for the system has been proposed in this work. Equations of motion have then been 
derived and transformed from the purely quantum mechanical description to a semi-classical picture 
using the method of coherent structures. We found that the dynamics of a MT can be explained by 
the cubic-quintic nonlinear Schrodinger equation (NLS) with a complex potential. A generic solution 
of the NLS equation in cylindrical geometry is a vortex filament [34-36] which can grow or shrink 
linearly in time as well as fluctuate temporally with some frequency. This behavior exhibits two 
distinct dynamical phases: (a) linear growth/shrinkage and (b) oscillation about a mean value, which 
are indeed main experimentally observed characteristics of the MT dynamics (Fig. 1 and 2). 

We have demonstrated here that the assembly process can be described starting from quantum me- 
chanical first principles applied to biochemical reactions. This can be subsequently transformed into a 
highly nonlinear semi-classical dynamics problem. The gross features of MT dynamics satisfy classical 
field equations in a coarse-grained picture. Individual chemical reactions involving the constituent 
molecules still retain their quantum character. The method of coherent structures allows for a simul- 
taneous classical representation of the field variables and a quantum approach to their fluctuations. 
Here, the overall MT structure (and their ensembles) can be viewed as a virtual classical object in 
(3+l)-D space-time. However, at the fundamental level of its constituent biomolecules, it is quantized 
as are chemical reactions involving its assembly or disassembly. In our approach the route taken is 
opposite since we started with individual tubulin quantum microstates to arrive at classical, nonlinear 
but classically coherent (and stable) macro-states of a microtubule. 

Finally, we note that dynamics of pattern formation can be also described by NLS equation in which 



al. [43] demonstrated that gravity can influence tubulin assembly reactions with MTs forming distinct 
bands at right angles to the orientation of the gravity field or, if spun, to the centrifugal force. Despite 
several studies [44-47] , the above experiments are yet to be fully explained theoretically. Our goal in 
future studies is to focus on the dynamics of pattern formation by MTs using the results presented in 
this paper. 
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